clear; clc;

m = 2000;
n = m;
k = floor(0.01*m);

X = randn(m, k)*randn(k, n);
maxX = max(abs(X(:)));
Y = 5*maxX*(sprand(m, n, 0.01) ~= 0);

D = X + Y;
O = D + 0.1*randn(size(X));

para.maxIter = 100;
para.tol = 1e-6;
para.objstep = 4;
para.tau = 1.01;
para.regType = 1;

lambdaMax = topksvd( O, 1, 5 );

%% start
gridLambda = lambdaMax*(1/2).^(2:10);
gridRMSE = zeros(size(gridLambda));
for g = 1:length(gridLambda)
    lambda = gridLambda(g);
    mu = lambda/sqrt(m);
    
    % Xa = APGRPCA( O, lambda, mu, para);
    Xa = FastProxRPCA( O, lambda, mu, 1e+8, 1e+8, para );
    
    gridRMSE(g) = sqrt(norm(Xa - X, 'fro')^2/numel(X));
    if(g > 1 && gridRMSE(g) > gridRMSE(g - 1))
        gridRMSE = gridRMSE(1:g);
        break;
    end
end

[~, lambda] = min(gridRMSE);
lambda = gridLambda(lambda);
mu = lambda/sqrt(m);

t = tic;
[ Xa, Ya, out{1} ] = APGRPCA( O, lambda, mu, para);
Time(1) = toc(t);
Da = Xa + Ya;
NMSE(1) = norm(Da - D, 'fro')/norm(D, 'fro');
Supp(1) = supportAcc( Ya, Y );
clear Xa Ya;

%% ---------------------------------------------------------------
for type = 1:3
para.regType = type;
para.decay = 0.6;
para.maxR = 2*k;
para.maxIter = 100;

if(type == 3)
    gridLambda = lambdaMax*(1/2).^(6:16);
else
    gridLambda = lambdaMax*(1/2).^(2:10);
end

gridRMSE = zeros(size(gridLambda));
for g = 1:length(gridLambda)
    lambda = gridLambda(g);
    
    switch (para.regType)
        case 1
            theta1 = lambda*2;
        case 2
            gridLambda(g) = 10*gridLambda(g);
            lambda = gridLambda(g);
            theta1 = sqrt(lambda);
        case 3
            theta1 = ceil(k/2);
    end

    Xa = FastProxRPCA( O, lambda, mu, theta1, 1e+8, para );
    
    gridRMSE(g) = sqrt(norm(Xa - X, 'fro')^2/numel(X));
    
    if(g > 1 && gridRMSE(g) >= gridRMSE(g - 1))
        gridRMSE = gridRMSE(1:g);
        break;
    end
end

[~, lambda] = min(gridRMSE);
lambda = gridLambda(lambda);

if(type ~= 3)
    gridTheta = 2.^(-1:2);
    gridRMSE = zeros(size(gridTheta));
    for g = 1:length(gridTheta)
        theta1 = gridTheta(g);
        switch (para.regType)
            case 1
                theta1 = theta1*lambda;
            case 2
                theta1 = theta1*sqrt(lambda);
            case 3
                theta1 = k;
        end

        Xa = FastProxRPCA( O, lambda, mu, theta1, 1e+8, para );
        gridRMSE(g) = sqrt(norm(Xa - X, 'fro')^2/numel(X));
    end

    [~, theta1] = min(gridRMSE);
    theta1 = gridTheta(theta1);
end

switch (para.regType)
    case 1
        theta1 = theta1*lambda;
    case 2
        para.maxR = k + 1;
        theta1 = theta1*sqrt(lambda);
    case 3
        lambda = lambda*2;
        theta1 = k;
end
para.maxIter = 500;

% ------------------------------------------------------------------------
t = tic;
[ Xa, Ya, out{2, type} ] = FastProxRPCA( O, lambda, mu, theta1, 1e+8, para );
Time(2, type) = toc(t);
Da = Xa + Ya;
NMSE(2, type) = norm(Da - D, 'fro')/norm(D, 'fro');
Supp(2, type) = supportAcc( Ya, Y );

% ------------------------------------------------------------------------
t = tic;
[ Xa, Ya, out{3, type} ] = GISTRPCA( O, lambda, mu, theta1, 1e+8, para );
Time(3, type) = toc(t);
Da = Xa + Ya;
NMSE(3, type) = norm(Da - D, 'fro')/norm(D, 'fro');
Supp(3, type) = supportAcc( Ya, Y );

clear Xa gridLambda gridRMSE gridTheta;

% ------------------------------------------------------------------------
% t = tic;
% [ Xa, Ya, out{4,type} ] = IRNNRPCA( O, lambda, mu, theta1, 1e+8, para );
% Time(4,type) = toc(t);
% Da = Xa + Ya;
% NMSE(4) = norm(Da - D, 'fro')/norm(D, 'fro');
% Supp(4,type) = supportAcc( Ya, Y );

end

clear X Y O D t n k Xa Ya g Da gridLambda gridRMSE gridTheta lambda lambdaMax;
clear maxX mu para theta1 type;

Supp = full(Supp);

save(strcat(num2str(m), '-rpca-1.mat'));
